1. Static vs dynamic data analysis

Static analysis is a common approach to formant analysis in phonetics research. It is often the case, for example, that we extract formant frequencies during a vocalic interval at a single point in time (e.g., segmental onset, temporal midpoint). This is a farily prevalent approach and it is useful when we try to understand general and abstract properties of segments.

While static analysis is a somewhat easy and computationally less expensive approach, this inevitably involves some degrees of data loss/compression. Recent studies have increasingly demonstrated that time-varying properties in speech signals provide rather fine-grained phonetic details that may be important to better understand e.g., synchronic and diachronic variation and change, and cross-linguistic differences/influences in the acquisition of phonetics/phonology. And this is where dynamic analysis methods will come in handy, in which we directly model non-linearity in the data without having to compress/losing much information. My motto in dynamic analysis is modelling what you see.

2. Top-down vs bottom-up approaches to statistical modelling

Another contrast that will be discussed in this workshop is the difference in the way we deal with hypotheses. On the one hand, we may be interested in seeing what the data show us when we do not have any prior knowledge. In other words, our study can be exploratory, and we sometimes prefer data-driven (or bottom-up) approach to identifying key variables of interest, which then feed hypothesis/research question generations.

On the other hand, statistical modelling can start with formulating hypotheses and research questions based on some existing ideas of key variables of interest, which could derive from e.g., previous research. In a way, this process involves a top-down decision making, and hypothesis testing involves asking whether the key variables of interest have a (statistically) significant effect in distinguishing two different populations.

The contrast between bottom-up (data-driven) and top-down (hypothesis-driven) approaches to data analysis influence the nature of research questions that can be asked using each method.

Summary

In the next two days, we will be mainly discussing four statistical methods corresponding to the contrasts above.

1. Principal Component Analysis (PCA) : static and bottom-up approach

2. Functional Principal Component Analysis (FPCA): dynamic and bottom-up approach

3. Linear Mixed-Effect Models (LMEs): static and top-down approach

4. Generalised Additive Mixed-Effect Models (GAMMs) : dynamic and top-down approach

Note that I do not wish to advocate one method over another: the main objective of this workshop is for us to experience that the same sets of data can be analysed differently, that data visualisation can look different, and that we can ask different research questions from the same sets of data.

In the next section, I will introduce a problem that we will be working on in the course of the workshop. Given the nature of my expertise, I only manage to contexualise the workshop on a somewhat specific research topic, but hopefully the workshop offers something exciting for everyone!

L1 Japanese speakers’ production of L2 English /l/ and /ɹ/

I contexualise this workshop broadly in the second language (L2) speech learning. More specifically, my research addresses the ‘’Japanese /r/-/l/ problem’’ (Flege, Aoyama & Bohn, 2021, p. 84). There is a long history of research on this topic, in which the first well-known, influential empirical studies could date back to 1970’s (e.g., Goto et al., 1971; Miyawaki et al., 1975). Up until today, lots of researchers have worked on this topic to explain the substantial difficulty that L1 Japanese speakers have in perceiving and producing L2 English liquids (e.g., Sheldon & Strange, 1982; Flege et al., 1995; Bradlow et al., 1997; 1999; Iverson et al., 2003; Aoyama et al., 2004; Saito & Munro, 2014; Saito & van Poeteren, 2018; Shinohara & Iverson, 2018; Aoyama et al., 2019; Aoyama et al., 2023).

Overarching question

There have been quite a few studies investigating L1 Japanese speakers’ production of L2 English liquids. To our surprise, however, there are not many studies in which the production data are analysed instrumentally using acoustic and/or articulatory methods. Previous instrumental acoustic studies would include: Flege et al. (1995), Saito and Munro (2014), Saito and van Poeteren (2018), Aoyama et al. (2019), Aoyama et al. (2023). This list is not exhaustive but does cover pretty much the main ones out there.

These studies commonly utilise static analysis in characterising the acoustic profiles of L1 Japanese speakers’ production of L2 English liquids. Thanks to these studies, we now understand the general acoustic properties involved in the production. These studies, however, have not managed to explain how and why L1 Japanese speakers end up producing non-target-like L2 English liquids.

So, in this workshop, we will explore ways to better explain the mechanism as to why L1 Japanese speakers have substantial difficulty in producing L2 English liquids in a target-like manner based on acoustic data. To rephrase this into a question:

What is the underlying mechanism of the L1 Japanese speakers’ difficulty in producing L2 English liquids in a target-like manner?

Workshop structure

In this two-day workshop, I hope to go through the whole procedure of data analysis from hypothesis generation to hypothesis testing. On Day 1, you will be faced with a new, unfamiliar data set, so we will explore the data structure and identify salient acoustic dimensions that are important in addressing the question. PCA and FPCA are strong tools at this stage for static and dynamic analyses. At the end of the day, ideally, we will be able to generate a hypothesis regarding between-group differences with specific acoustic parameters. On Day 2, we will model and test the hypothesis generated on Day 1 using LMEs and GAMMs.

Workshop materials

In this workshop, we will be using the materials from one of my previous publications. Data and codes for the analysis are publicly available on the Open Science Framework (OSF) repository, and we will retrieve the data sets from there.

If interested, the article can be accessed here (open access):

Nagamine, T. (2024). Formant dynamics in second language speech: Japanese speakers’ production of English liquids. The Journal of the Acoustical Society of America, 155(1), 479–495. https://doi.org/10.1121/10.0024351

Here is the URL for the online supplementary materials deposited on the following OSF repository:

https://osf.io/2phx5/

Mean and median: the distribution “centre”

The most common parameters in describing a given set of data are mean and median. Assuming the normal distribution, the most common type of data distribution in statistics, both mean and median represent something like the midpoint. But they differ in a small detail:

• Mean: the sum of all values divided by the number of values. Sensitive to extreme values.

• Median: the middle value when lining up all the data points from the smallest to the greatest. More robust to extreme values than the mean.

The figures below illustrate mean and median values in (1) normal and (2) right skewed distributions. You can see that mean and median values match on the normal distribution (left plot), whereas the mean score is ‘attracted’ towards the extreme values towards the left (right plot).

Standard deviation: the distribution “spread”

Along with the mean and median values, we often encounter standard deviation (SD). SD captures the degree of ‘’spread’’ of a given distribution, or the distance from the mean in the distribution.

SD tells us the degree of data coverage. For example, one of the distributions shown below has the mean of 50 and the SD of 30 (in the skyblue colour). This indicates that 68% of the data points are covered between -1SD and +1SD: that is, between 50 - 30 = 20 and 50 + 30 = 80.

Similarly, 95% of the data points are included within ±2SDs. In the skyblue distribution below, 95% of the data are covered between 50 - 30 \(\times\) 2 = -10 and 50 + 30 \(\times\) 2 = 110. Usually, the larger SD is, the more spread the distribution looks.

Linear vs non-linear data

Let’s take a quick look into what it means by ‘linear’ and ‘non-linear’ data. Both are the terms used to explain the relationship between two variables. A linear relationship can be expressed as a straight line, whereas a non-linear relationship requires a curve (or a contour).

The figure below exemplifies a linear and a non-linear relationship between X and Y. Generally, the plot with dots like this is called ‘scatter plots’. Here, the black dots show a linear relationship, where a certain amount of increase in X corresponds to another certain amount of increase in Y, and this X-Y correspondence is consistent throughout the range of X/Y values. This linear relationship can be expressed using a straight (black) line.

The yellow dots, on the other hand, show a non-linear relationship between X and Y variables. Although an increase in X corresponds to another increase in Y like in the linear data, the relationship is not consistent, especially towards the end of the x-axis. This relationship cannot be captured with a single line (compare the solid and dotted lines in yellow).

Data visualisation: scatter, box and violin plots

Finally, let’s take a quick look at data visualisation methods. I consider that data visualisation is the pivot in data analysis, and a great data analysis is facilitated by effecive data visualiastion. It is also a means to communicate the research findings with the audience.

Effective data visualisation often combines multiple types of plots. In my usual practice, I often combine three types of plots: scatter, box and violin plots.

1. Scatter plot: In the figure below, the dots in skyblue colour represents individual data points. The mean value of the data is expressed as ▲.

2. Box plot: Superimposed on them is the box plot. The box plot consists of a rectangle box, a solid horizontal line in the middle of the box, and ‘whiskers’ at the top and bottom of the box. Each of them represents key parameters of a distribution:

• The solid horizontal line represents median.

• Each end of the box represents first and third quartiles. Broadly speaking, you could understand that the box covers the middle 50% of the data, which is divided into half by the solid black line. The lower end of the box is the first quartile (Q1) and the upper end is the third quartile (Q3). Right in the middle is the median, which is the second quartile (Q2).

• The length of the whiskers below and above the box is determined such that:

  • Lower end: Q1 - 1.5 \(\times\) IQR

  • Upper end: Q3 + 1.5 \(\times\) IQR

3. Violin plot: Finally, the violin plot shows the density of the data distribution.

Installing/loading packages

Let’s first install and load R packages that we are using in the static analysis section. The installation commands have been commented out, but please uncomment them and install the packages if you do not have them on your machine yet. We won’t need any specialised package for PCA as it is included as a base function.

# installing
# install.packages("tidyverse")
# install.packages("lme4")
# install.packages("lme4Test")
# install.packages("emmeans")

# importing
library(tidyverse)
library(lme4)
library(lmerTest)
library(emmeans)

Importing data set

Let’s import the data set. We are using the data set openly available on the Open Science Framework (OSF) repository.

# import the csv file "initial.liquid.static.csv" from the "data" directory and save it as "df_mid"
df_mid <- readr::read_csv("data/initial.liquid.static.csv")

Check data

First, we will inspect a summary of the data according to the columns that we have so far. Although not the most efficient, my approach is to combine dplyr::group_by() and dplyr::summarise() functions to summarise the data.

Do not forget dplyr::ungroup() at the end of each operation, too!

# Let's check the number of "approximant" tokens
df_mid |> 
  dplyr::group_by(IsApproximant) |> 
  dplyr::summarise() |> 
  dplyr::ungroup()

## # A tibble: 1 × 1
##   IsApproximant
##   <chr>        
## 1 yes

# Let's check the number of tokens of good recording quality
df_mid |> 
  dplyr::group_by(IsAcoustic) |> 
  dplyr::summarise() |> 
  dplyr::ungroup()

## # A tibble: 1 × 1
##   IsAcoustic
##   <chr>     
## 1 yes

# How about 'omit'?
df_mid |> 
  dplyr::group_by(omit) |> 
  dplyr::summarise() |> 
  dplyr::ungroup()

## # A tibble: 1 × 1
##    omit
##   <dbl>
## 1     0

Omitting irrelavent columns

The codes above show that the data set only contains liquid tokens that are classified as an approximant, deemed as being of good recording quality and thus included in the acoustic analysis. This means that we can safely remove these columns. Let’s do this via dplyr::select() function. After that, we check the column names again using colnames(), showing that the three columns have been safely removed from the data frame.

# Remove columns 
df_mid <- df_mid |> 
  dplyr::select(-c(IsApproximant, IsAcoustic, omit))

# Check column names again
# ...

Renaming variables

Later on in this workshop, we will be modelling the formant frequencies as a function of vowel context. In terms of the current data set, this corresponds to the next_sound column, but let’s rename it into vowel to make it more intuitive.

Also, we will also convert the ARPABET notations into the IPA symbols.

df_mid <- df_mid |> 
  dplyr::mutate(
    vowel = case_when(
      next_sound == "AE1" ~ "/æ/",
      next_sound == "IY1" ~ "/i/",
      next_sound == "UW1" ~ "/u/",
    )
  )

Checking the number of participants, tokens…

Let’s also obtain some descriptive statistics here.

# number of participants
df_mid |> 
  dplyr::group_by(language) |> 
  dplyr::summarise(n = n_distinct(speaker)) |> 
  dplyr::ungroup()

## # A tibble: 2 × 2
##   language     n
##   <chr>    <int>
## 1 English     14
## 2 Japanese    31

# number of tokens per segment
df_mid |> 
  dplyr::group_by(segment) |> 
  dplyr::summarise(n = n()) |> 
  dplyr::ungroup()

## # A tibble: 2 × 2
##   segment     n
##   <chr>   <int>
## 1 L        1114
## 2 R        1192

# Check data further
# ...

Preparation

Before visualising the data, we need to make sure that the formant frequencies are comparable across speakers. We know that we cannot directly compare raw formant frequencies because they are heavily influenced by the vocal tract characteristics depending on the speaker’s gender, sex, height etc.

To address this, let’s normalise the formant frequencies across speakers. There are various normalisation methods, but today we are just using z-score normalisation, in which the mean value for each speaker is expressed as zero and the standard deviation is scaled to one. This can be achieved via scale() function included in base R.

Let’s create f1z, f2z and f3z columns using the dplyr::mutate() function. When doing this, make sure that you do the scaling for each speaker by specifying dplyr::group_by(speaker).

df_mid <- df_mid |> 
  dplyr::group_by(speaker) |> # tell R to do the following iteration per speaker
  dplyr::mutate(
    f1z = as.numeric(scale(f1)), # scale f1 into z-score
    f2z = as.numeric(scale(f2)), # scale f2 into z-score
    f3z = as.numeric(scale(f3)) # scale f3 into z-score
  ) |> 
  dplyr::ungroup() # don't forget ungrouping

Let’s check the mean and SD for both raw and normalised formant values: just see F1 for now. Note that the mean z-scores do not seem to look zero, but this is because computers are not very good at dealing with very small numbers (e.g., decimals) and some fluctuations occur in computing the values.

# check mean and sd of raw/scaled F1 values for each speaker
df_mid |> 
  dplyr::group_by(speaker) |>
  dplyr::summarise(
    f1_mean = mean(f1),
    f1_sd = sd(f1),
    f1z_mean = mean(f1z),
    f1z_sd = sd(f1z)
  ) |> 
  dplyr::ungroup() 

## # A tibble: 45 × 5
##    speaker f1_mean f1_sd  f1z_mean f1z_sd
##    <chr>     <dbl> <dbl>     <dbl>  <dbl>
##  1 2d57ke     372.  74.7  3.13e-16      1
##  2 2drb3c     430.  50.6  3.04e-16      1
##  3 2zy9tf     358.  73.8 -3.20e-16      1
##  4 3bcpyh     374.  41.0  1.68e-16      1
##  5 3hsubn     390.  82.6 -2.39e-16      1
##  6 3pzrts     335.  27.7  6.52e-16      1
##  7 4ps8zx     392.  68.7 -3.06e-16      1
##  8 54i2ks     348.  65.9 -3.35e-17      1
##  9 5jzj2h     366.  76.8 -1.22e-16      1
## 10 5upwe3     343.  42.1  1.05e-16      1
## # ℹ 35 more rows

# check mean and sd of raw/scaled F2 values for each speaker group
# ...

# check mean and sd of raw/scaled F3 values for each speaker group
# ...

Visualisation: Spectral characteristics

We can now compare F1 values across participant groups for each segment (/l/ and /ɹ/). For clearer visualisation, the following codes define the range of y-axis between -4 and 4, and some extreme values have been omitted from the plot.

In the code below, I visualise the f1z values using the scatter, box and violin plots. The line with ggplot() defines x- and y-axis, and the following three lines with geom_jitter(), geom_violin() and geom_boxplot() superimpose the plot on top of each other and define aesthetic. facet_grid() divides the plot into multiple levels of a variable to understand differences across conditions.

# F1
df_mid |> 
  # define x- and y-axis
  ggplot(aes(x = language, y = f1z)) +
  # scatter plot
  geom_jitter(aes(colour = language), width = 0.3, alpha = 0.4) +
  # superimposing violin plot on scatter plot
  geom_violin(alpha = 0.4) +
  # superimposing boxplot on scatter/violin plots
  geom_boxplot(width = 0.4, alpha = 0.4) +
  # drawing a horizontal line at y = 0 to facilitate interpretation
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.5) +
  # define colour scheme
  scale_colour_manual(values = cbPalette) + 
  # define limits on the y-axis
  ylim(c(-4, 4)) +
  # split the plot according to the ``segment`` and ``vowel`` variables
  facet_grid(segment ~ vowel) +
  # define axis and main labels
  labs(x = "L1", y = "F1 (z-normalised)", title = "F1 frequency") +
  # adjust the label on the y-axis to make it look better
  theme(strip.text.y = element_text(angle = 0))

# F2
# df_mid |> ...
 
# F3
# df_mid |> ...

Exploring visualisation further

# try your visualisation
# ...

Correlation

Let’s also take a slightly different approach. The objective of PCA is to draw straight lines along the dimension that shows the greatest variance. To inspect the variance in the data, we’ll plot the relationships between the three spectral dimensions.

# F1 and F2
df_mid |> 
  ggplot(aes(x = f1z, y = f2z)) +
  geom_point(aes(colour = segment), alpha = 0.4) +
  scale_colour_manual(values = cbPalette) +
  facet_wrap(~ language)

# F2 and F3
df_mid |> 
  ggplot(aes(x = f2z, y = f3z)) +
  geom_point(aes(colour = segment), alpha = 0.4) +
  scale_colour_manual(values = cbPalette) +
  facet_wrap(~ language)

# F1 and F3
df_mid |> 
  ggplot(aes(x = f1z, y = f3z)) +
  geom_point(aes(colour = segment), alpha = 0.4) +
  scale_colour_manual(values = cbPalette) +
  facet_wrap(~ language)

As expected, English /l/ and /ɹ/ can be reliably distinguished along the F3 dimension, such that English /ɹ/ shows lower F3 than English /l/.

Temporal characteristics

OK, then how about duration? Let’s visualise them first, and then check how it correlates with each of the spectral parameters. We’ll use the original df_mid data frame again but convert the duration into millisecond.

# convert duration (s) into (ms)
df_mid <- df_mid |> 
  dplyr::group_by(speaker) |> 
  dplyr::mutate(
    duration_ms = duration * 1000,
    duration_ms_z = scale(duration_ms)
  ) |> 
  dplyr::ungroup()

# plot
df_mid |> 
  ggplot(aes(x = language, y = duration_ms_z)) +
  geom_jitter(aes(colour = language), alpha = 0.5) +
  geom_violin(alpha = 0.4) +
  geom_boxplot(width = 0.3, alpha = 0.4) +
  facet_grid(~ segment) +
  scale_colour_manual(values = cbPalette)

# with F1 
df_mid |> 
  ggplot(aes(x = f1z, y = duration_ms_z)) +
  geom_point(aes(colour = language), alpha = 0.4) +
  scale_colour_manual(values = cbPalette) +
  facet_grid(segment ~ language) +
  theme(strip.text.y = element_text(angle = 0))

# with F2
df_mid |> 
  ggplot(aes(x = f2z, y = duration_ms_z)) +
  geom_point(aes(colour = language), alpha = 0.4) +
  scale_colour_manual(values = cbPalette) +
  facet_grid(segment ~ language) +
  theme(strip.text.y = element_text(angle = 0))

# with F3
df_mid |> 
  ggplot(aes(x = f3z, y = duration_ms_z)) +
  geom_point(aes(colour = language), alpha = 0.4) +
  scale_colour_manual(values = cbPalette) +
  facet_grid(segment ~ language) +
  theme(strip.text.y = element_text(angle = 0))

PCA summary

Let’s also look into the inside of pca_mid.

# show the PCA results
pca_mid

## Standard deviations (1, .., p=4):
## [1] 1.1419985 1.0410699 0.9836312 0.7537420
## 
## Rotation (n x k) = (4 x 4):
##                      PC1          PC2        PC3        PC4
## f1z           0.02567994  0.554175433  0.7973415  0.2376484
## f2z           0.59519981  0.375814985 -0.4453013  0.5533598
## f3z           0.73795645 -0.003335683  0.1729687 -0.6522967
## duration_ms_z 0.31701420 -0.742725764  0.3688295  0.4602447

The top line shows standard deviation, which we just saw earlier, so let’s skip this for now.

At the bottom, we have loadings matrix. This shows eigenvectors: i.e., the amount and (2) the direction of the contribution that the original data dimensions make to each PC. Let’s disentangle this one by one:

• PC1: Contributions of f1z, f2z, f3z and duration are all in the same direction (i.e., positive). But the amount of contribution is fairly small for f1z and duration compared to f2z and f3z.

• PC2: In contrast to PC1, f1z contributes to PC2 to a greater extent than f2z and f3z. duration also seems to make a big contribution here but in the opposite direction (i.e., negative).

• PC3: The contribution of duration seems quite big here, too.

• PC4: f2z and f3z both seem to contribute to this dimension, although they are in antagonistic direction.

So, all in all, PCA might suggest:

• PC1 captures a covariation between f2z and f3z, which may correspond to the spectral difference between /l/ and /ɹ/.

• PC2 mainly captures a covariation of duration and f1z but in opposite directions.

• PC3 shows variation in f1z.

• PC4 again seems to capture f2z and f3z but in opposite directions.

The interpretation may be better facilitated by data visualisation. Let’s try two main approaches to data visualiastion of PCA.

Data visualisation 1: Scree plot

The amount of variance is useful when determining which PC to retain for analysis. As a rule of thumb, Bayeen (2008) recommends that we retain PCs that explains greater than 5% of variance in the data.

A useful visualisation of the proportion of variance is scree plot. There is a default function screeplot() that lets you create a scree plot quite easily:

# scree plot: bar chart
screeplot(pca_mid)

# scree plot: line plot
screeplot(pca_mid, type = "line")

You can also create a scree plot using ggplot() by manually calculating the proportion of variance based on the standard deviation (see above for the formula). This is useful when you need more flexibility – e.g., to show the 5% thereshold suggested by Bayeen (2008).

## obtaining proportion of variance 
pca_var_exp <- pca_mid$sdev^2 / sum(eigenvalues)

# making var_explained as a tibble 
pca_var_exp <- tibble::as_tibble(pca_var_exp)

# adding column name
pca_var_exp <- pca_var_exp |>  
  tibble::as_tibble() |>  
  dplyr::mutate(
    PC = row_number()
  )

# create a plot
scree_plot <- pca_var_exp |> 
  ggplot(aes(x = PC, y = value)) +
  geom_line() +
  geom_text(aes(label = round(value, digits = 3)*100), nudge_x = 0, nudge_y = -0.02) +
  geom_label(aes(label = PC), label.padding = unit(0.40, "lines")) +
  geom_hline(yintercept = 0.05, linetype = 'dotted') +
  scale_x_continuous(breaks = c(1, 2, 3)) +
  labs(x = "Principal Component", y = "Variance Explained", title = "Proportion of Variance explained by each PC")

# showing plot
scree_plot

Since there are only three PCs identified from the data, it doesn’t make much sense to visualise them. But this will be useful when prcomp() identifies more than a few PCs.

creating biplot

The interpretation of loading matrix was more or less straightforward in this case given that only three parameters were involved. It is easy to see, however, that this can get quite complicated when the original data have a greater number of dimensions.

Biplots help us better understand the loadings and PC scores for each observation. For example, the default function biplot() plots all observations along the PC1 and PC2 dimensions with eigenvectors shown in red arrows (i.e., the direction and the amount of weighting: also called loadings).

# default function for biplot
biplot(pca_mid)

Detailed biplots

Similarly to the scree plot, we could also create biplots manually using ggplot(). This will help us better understand what each PC could stand for in the current data set.

# Extract PCA scores (principal component scores for each observation)
pca_scores <- as.data.frame(pca_mid$x)

# Extract PCA loadings (the contribution of each original variable to the PCs)
pca_loadings <- as.data.frame(pca_mid$rotation)

# Combine scores into one data frame for plotting
scores_df <- cbind(pca_scores, df_mid[, c("segment", "language", "vowel")])

# Reshape the loadings for visualization
loadings_df <- pca_loadings |> 
  mutate(variable = rownames(pca_loadings))

# Plotting the PCA biplot with facets by language
ggplot() +
  geom_point(data = scores_df, aes(x = PC1, y = PC2, color = segment, shape = segment), size = 2, alpha = 0.5) + # Plot the PCA scores (observations)
  geom_segment(data = loadings_df, aes(x = 0, y = 0, xend = PC1 * 5, yend = PC2 * 5), arrow = arrow(type = "closed", length = unit(0.1, "inches")), color = "red") + # Plot the loadings (variables)
  geom_text(data = loadings_df, aes(x = PC1 * 5, y = PC2 * 5, label = variable), size = 4, color = "red") + # Add labels for the loadings
  # stat_ellipse(data = scores_df, aes(x = PC1, y = PC2, color = vowel, fill = vowel), geom = "polygon", alpha = 0.2, level = 0.95) + 
  facet_grid(language ~ vowel) +
  labs(title = "PCA biplot by segment and vowel") +
  scale_color_manual(values = cbPalette) + 
  scale_fill_manual(values = cbPalette) + 
  theme(legend.position = "bottom",
        strip.text.y = element_text(angle = 360))

Extracting PC scores

PC scores are stored as x in the PCA output. Let’s use head() function to only display the first six rows as otherwise the list would be too big.

# displaying PC scores
head(pca_mid$x)

##             PC1         PC2         PC3        PC4
## [1,]  0.1030871 -0.46456046 1.842069714 -0.3630188
## [2,]  0.3992679 -0.34912107 1.606508257 -0.3671703
## [3,] -0.2207061  0.81275522 1.210015271 -1.0763739
## [4,]  0.1631096 -0.02578864 0.004715655 -1.2024286
## [5,]  0.2346244 -0.63551312 1.461971228 -0.9505229
## [6,]  0.3949973 -0.74218923 2.054526744 -0.1311160

Assuming that we haven’t made any changes to the order of rows, we can combine the PC scores with the existing data set. We have actually already done this when visualising the data, so I simply copy and paste the codes here:

# extract PCA scores (principal component scores for each observation)
pca_scores <- as.data.frame(pca_mid$x)

# check the number of rows
## pc score
nrow(pca_scores)

## [1] 2306

## main data set
nrow(df_mid)

## [1] 2306

Both pc_score and df_mid contain the same number of rows, so we can merge them into one data frame.

# combine scores into one data frame for plotting
scores_df <- cbind(pca_scores, df_mid[, c("segment", "language", "vowel", "speaker", "gender", "word")])

Check PC scores

Before visualisation, let’s take a look at descriptive statistics of PC scores. We’ll focus on the first two PCs only here.

scores_df |> 
  dplyr::group_by(language, segment, vowel) |> 
  dplyr::summarise(
    PC1_mean = mean(PC1),
    PC1_sd = sd(PC1),
    PC2_mean = mean(PC2),
    PC2_sd = sd(PC2)
  ) |> 
  dplyr::ungroup()

## `summarise()` has grouped output by 'language', 'segment'. You can override
## using the `.groups` argument.

## # A tibble: 12 × 7
##    language segment vowel PC1_mean PC1_sd PC2_mean PC2_sd
##    <chr>    <chr>   <chr>    <dbl>  <dbl>    <dbl>  <dbl>
##  1 English  L       /i/      0.917  0.785 -0.0570   0.923
##  2 English  L       /u/      0.965  0.762 -0.629    0.989
##  3 English  L       /æ/      0.797  0.593  0.104    1.08 
##  4 English  R       /i/     -0.632  0.768  0.0538   0.940
##  5 English  R       /u/     -0.808  0.685 -0.177    0.978
##  6 English  R       /æ/     -1.06   0.769  0.446    1.33 
##  7 Japanese L       /i/      0.898  0.753 -0.00896  0.927
##  8 Japanese L       /u/      0.278  0.956 -0.617    0.953
##  9 Japanese L       /æ/      0.585  0.890  0.146    0.995
## 10 Japanese R       /i/     -0.163  1.03   0.221    0.905
## 11 Japanese R       /u/     -0.901  0.907 -0.392    0.885
## 12 Japanese R       /æ/     -0.721  1.08   0.218    1.03

We get lots of numbers here – please do feel free to take a moment to digest these. But we could also visualise them to better facilitate the interpretation.

Visualising PC1

Let’s visualise PC scores to see if our assumptions about each PC is correct. We think that PC1 captures covariation of f2z and f3z, and this should correspond to the contrast between English /l/ (lower F2, higher F3) and /ɹ/ (higher F2, lower F3).

## PC1 - by liquid consonant
scores_df |> 
  ggplot(aes(x = segment, y = PC1, colour = segment)) +
  geom_jitter(width = 0.3) +
  geom_violin(alpha = 0.3) +
  geom_boxplot(width = 0.3, colour = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.3) +
  scale_colour_manual(values = cbPalette) +
  facet_grid(language ~ vowel) +
  theme(
    strip.text.y = element_text(angle = 360)
  )

This is more or less true, good! We could also make a between-group comparison.

## PC1 - by group
scores_df |> 
  ggplot(aes(x = language, y = PC1, colour = language)) +
  geom_jitter(width = 0.3) +
  geom_violin(alpha = 0.3) +
  geom_boxplot(width = 0.3, colour = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.3) +
  scale_colour_manual(values = cbPalette) +
  facet_grid(segment ~ vowel) +
  theme(
    strip.text.y = element_text(angle = 360)
  )

We can see some patterns here: e.g., both L1 English and L1 Japanese speakers are similar for /l/, L1 Japanese speakers produce smaller PC1 scores in the /u/ context than L1 English speakers. Some between-group differences can be found for English /ɹ/ as well.

Visualising PC2

## PC2
### scores_df |> ...

Statistical analysis

As I said earlier, as a somewhat common practice, PC scores can serve as input to further statistical analysis. So far, we have found that PC1 captures the contrast between English /l/ and /ɹ/; although this is sort of obvious, the most important message here is that the data tells us about it without any apriori knowledge.

In the data visualisation above, it was somewhat difficult to identify which of these variables: language, vowel and segment, would have a significant effect. At this point, we could construct a linear mixed-effect model to test whether these variables have a statistically significant effect on the PC scores.

I don’t show the full results here because it’s too long, but more to come tomorrow when we come back to talk about LMEs!

# converting variables into factor
scores_df <- scores_df |> 
  dplyr::mutate(
    language = as.factor(language),
    vowel = as.factor(vowel),
    segment = as.factor(segment),
    speaker = as.factor(speaker),
    word = as.factor(word)
  )

# run a full model
m1 <- lme4::lmer(PC1 ~ language + segment + vowel + language:segment + language:vowel + (1|speaker) + (1|word), data = scores_df, REML = FALSE)

# model summary
summary(m1)

# nested model 1 -- testing the language-segment interaction
m2 <- lme4::lmer(PC1 ~ language + segment + vowel + language:vowel + (1|speaker) + (1|word), data = scores_df, REML = FALSE)

# model comparison
anova(m1, m2, test = "Chisq")

# nested model 2 -- testing the language-vowel interaction
m3 <- lme4::lmer(PC1 ~ language + segment + vowel + language:segment + (1|speaker) + (1|word), data = scores_df, REML = FALSE)

# model comparison
anova(m1, m3, test = "Chisq")

# post-hoc analysis
## between-group effects of each language
emmeans::emmeans(m1, pairwise ~ language | segment:vowel, adjust = "tukey")

Saving the model

It will be interesting to fit a linear mixed-effect model on PC scores tomorrow. So let’s save the model scores_df so that we can retrieve it easily when we need it.

save(scores_df, file = "data/pca_scores_df.rda")